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ABSTRACT 

We combine dynamical and non-equilibrium chemical modeling of evolving prestellar molecular 
cloud cores and investigate the evolution of molecular abundances in the contracting core. We model 
both magnetic cores, with varying degrees of initial magnetic support, and non-magnetic cores, with 
varying collapse delay times. We explore, through a parameter study, the competing effects of various 
model parameters in the evolving molecular abundances, including the elemental C/O ratio, the 
temperature, and the cosmic-ray ionization rate. We find that different models show their largest 
quantitative differences at the center of the core, whereas the outer layers, which evolve slower, have 
abundances which are severely degenerate among different dynamical models. There is a large range 
of possible abundance values for different models at a fixed evolutionary stage (central density), which 
demonstrates the large potential of chemical differentiation in prestellar cores. However, degeneracies 
among different models, compounded with uncertainties induced by other model parameters, make it 
difficult to discriminate among dynamical models. To address these difficulties, we identify abundance 
ratios between particular molecules, the measurement of which would have maximal potential for 
discrimination among the different models examined here. In particular, we find that the ratios 
between NH3 and CO; NH2 and CO; NH3 and HCO + are sensitive to the evolutionary timescale, 
and that the ratio between HCN and OH is sensitive to the C/O ratio. Finally, we demonstrate that 
measurements of the central deviation (central depletion or enhancement) of abundances of certain 
molecules are good indicators of the dynamics of the core. 

Subject headings: ISM: molecules - ISM: clouds - ISM: dust - magnetic fields - MHD - stars: formation 
- ISM: abundances 
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1. INTRODUCTION 

The earliest phases of star formation have been the 
subject of intense observational and theoretical investi- 
gation over the past four decades. Nevertheless, together 
with the very late stages of star formation (formation of 
planetary systems), the very early processes of molec- 
ular cloud fragmentation into high-density "cores" and 
the subsequent early contraction of these cores before 
protostars form in their centers remain the two least con- 
strained and most debated aspects of star formation the- 
ory (McKee and Ostriker 2007). The reason for this lack 
of strong constraints is that it is difficult to extract phys- 
ical conditions in star forming regions from observable 
quantities, such as dust column or molecular line maps; 
each observable is typically affected by multiple factors 
which are nontrivial to deconvolve. As a result, theo- 
retical models of cloud fragmentation and core evolution 
span a large range of possibilities. 

A diverse array of mutually interacting microphysical 
processes influence the evolution of molecular clouds to 
differing extents: gravity, turbulence, magnetic fields, 
UV and cosmic ray (CR) ionization, and a complex net- 
work of chemical reactions. Without strong observa- 
tional constraints on the relative importance of these 
processes, competing star formation theories have arisen 
based on different assumptions of the physics driving the 
formation of prestellar cores and their collapse to form 
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protostars. 

The path to the appearance of a protostar depends 
sensitively on the initial conditions and the relative im- 
portance of the different physical processes driving the 
formation and contraction of prestellar molecular cloud 
cores. The most important distinction among different 
calculations of prestellar collapse pertains to the nature 
of any forces opposing gravity during the prestellar col- 
lapse. In this context, we can identify two broad classes 
of models: non-magnetic and magnetic. 

Non-magnetic collapsing models include early numeri- 
cal and semi- analytic solutions of the Euler equations of 
hydrodynamic collapse (e.g., Larson 1969; Penston 1969; 
Yorke & Kriigel 1977; Shu 1977) 3 , as well as numerical 
studies of the gravo-turbulent formation and evolution of 
bound condensations formed by converging flows in tur- 
bulent model clouds (see Mac Low & Klessen 2004 and 
Scalo & Elmegreen 2004 for recent reviews). The evolu- 
tion of cores formed by turbulence-driven fragmentation 
is compatible with models of pure hydrodynamical col- 
lapse (e.g., Gong & Ostriker 2009). Turbulence-initiated 
core formation models can be generally separated into 
models of rapid global collapse, where the entire molecu- 
lar cloud is viewed as an evanescent structure (e.g., Hart- 
mann et al. 2001; Heitsch & Hartmann 2008), and mod- 
els with some initial decaying turbulent support that ex- 
tends the lifetime of clouds (e.g., Krumholz et al. 2006; 

3 Note however that Shu's singular isothermal sphere models 
refer to protostellar cores in which a central object exists, since 
the initial conditions for the inside-out collapse are singular at the 
center. Thus, this class of models cannot be used to study the early 
prestellar phase of the core evolution. 
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Tan et al. 2006). 

Magnetic models include magnetohydrodynamical cal- 
culations of the magnetically mediated contraction of 
prestellar cores, which is typically slower than the pure 
hydrodynamical collapse (with lower infall velocities and 
longer evolution timescales, e.g. Fiedler & Mouschovias 
1993; Li & Shu 1996; Tassis & Mouschovias 2007), and 
result in a characteristic hourglass-like configuration of 
the magnetic field threading the core (see Mouschovias & 
Ciolck 1999 for a review). In magnetic models, the life- 
time of the parent molecular cloud is controlled by the 
initial value of the mass-to-magnetic-flux ratio. If this 
ratio is smaller than a critical value (the cloud is magnet- 
ically subcritical), then the cloud as a whole is long-lived 
as the magnetic field can support it against gravity. In 
this case, magnetically supercritical cores form through 
the process of ambipolar diffusion and collapse to form 
protostars, while the cloud envelope remains supported 
by the field. If the mass to magnetic flux ratio of the 
parent cloud is larger than the critical value, then the 
cloud as a whole is unstable and short-lived, while the 
formation and collapse of cores also proceeds in shorter 
timescales. 

The coupling between chemistry and dynamics is an 
essential step in connecting the dynamics of prestellar 
cores to observations, and has long been recognized as 
critical in extracting information from molecular line ob- 
servations. The problem has been pursued for the past 
decade in a series of ground-breaking studies. Bergin 
& Langer (1997) identified, through studies of static 
media, "early" and "late" molecules (molecules which 
peak in abundance either in the initial or in the later 
stages of the evolution of a molecular cloud fragment). 
Aikawa et al. (2001) compared observations of L1544 
with coupled chemical/dynamical models. For the dy- 
namics, they used the analytical similarity solution of 
Larson (1969) and Penston (1969), as well as a "slowed 
down" version of it (using an overall slow-down factor) 
to approximate the slower magnetic collapse, and found 
that their fast models are better fits to L1544. Li et al. 
(2002) and Shematovich (2003) also modeled L1544, us- 
ing a non-magnetic model and the Li (1999) magnetic 
model for the dynamics, which treats the magnetic field 
in an approximate way including only magnetic pressure 
forces and assuming a spherical geometry. They also in- 
cluded grain-surface reactions in their calculations, and 
they found the magnetic model to be a better fit to the 
observations. They showed that the magnetic collapse 
is accelerating, and therefore a uniform slow-down of 
the non-magnetic collapse is not a good approximation, 
because in magnetic collapse the amount of time spent 
at different densities changes with respect to the non- 
magnetic case. Aikawa et al. (2003) updated the Aikawa 
et al. (2001) calculation to include grain-surface reac- 
tions, and they found that moderate retardation factors 
(factor of 3) are not in as strong disagreement with the 
data as in Aikawa et al. (2001), however they still found 
that the faster models are preferred. Li et al. (2002) 
however noted that their preferred Larson-Penston solu- 
tion features very large infall velocities that cannot be 
reconciled with observations of L1544. Lee et al. (2003) 
modeled the chemistry of various prestellar cores to fit 
Bonnor-Ebert and Plummer-like radial density profiles 
using dust continuum and line emission. Lee et al. (2004) 



looked at the chemical evolution of prestellar cores aim- 
ing to produce line profiles for prestellar and protostcllar 
cores. For modeling the dynamical evolution of prestel- 
lar cores they used a sequence of Bonnor-Ebert profiles, 
however the assignment of times and velocities to each 
profile was not treated self-consistently. For the proto- 
stellar cores they used the inside-out collapse model of 
Shu (1977). Flower et al. (2005), using a homologous col- 
lapse treatment of prestellar cores, emphasized the effect 
of grain growth on the chemistry. Aikawa et al. (2005) 
revisited the prestellar core problem, but instead of using 
an analytic solution for the dynamics, they used post- 
processing of a numerical simulation of the collapse of 
Bonnor-Ebert spheres (one very close to equilibrium and 
one very unstable), and found significant differences in 
the system chemistry depending on the initial conditions. 
This result was confirmed by Keto & Caselli (2008) who 
coupled one-dimensional dynamical simulations with a 
simple CO chemistry network and radiative transfer. In 
a follow-up study, Keto & Caselli (2010), starting from 
a Bonnor-Ebert sphere of central density 2 x 10 cm~ 3 
in unstable equilibrium, evolving quasi-statically, found 
that at a density of 2 x 10 7 cm -3 the abundance of CO is 
much lower than observations in LI 544, and argued for 
a higher CO desorption rate from grains. 

On certain issues a consensus is forming among these 
studies (see, e.g., Bergin & Tafalla 2007 for a recent re- 
view). One such issue is the sensitivity of the chemical 
state of the contracting core to the duration of the col- 
lapse. Characteristic examples are the carbon-bearing 
molecules, which tend to become depleted on the ice 
mantles of dust grains and are therefore most prominent 
in the early stages of evolution of a core, in contrast to 
nitrogen-bearing molecules which are not depleted and, 
as a result, their abundance increases in later stages. 

However, core age is not the only contributing fac- 
tor in determining molecular abundances. Observations 
have established the existence of chemical differentiation 
in prestellar cores of otherwise similar properties, which 
however show strong variations in the abundances of a 
number of chemical species. Characteristic examples of 
chemical differentiation can be seen in observations of 
prestellar cores L1544, L1498, L1517B which show CO 
depletion (Tafalla et al. 2002, 2004), while L1521E (a 
core of similar central density with L1517B) shows no 
CO depletion (Tafalla & Santiago 2004). This differenti- 
ation may be a result of the dynamics of the collapse as 
suggested by, e.g., Shematovich et al. (2003), and thus 
reflect the initial conditions in the parent star forming 
region. 

In this paper, we self-consistently follow the chemistry 
and dynamics of contracting prestellar cores. We do so 
by adding to the system of equations of Eulerian hy- 
drodynamics, or MHD (as appropriate), one continuity 
equation with sources and sinks given by the chemical 
reactions network for each species followed. In this way, 
the dynamics and non-equilibrium chemistry evolve si- 
multaneously. This is particularly important in the case 
of MHD models, as abundances of charged species have a 
feedback effect on the dynamics, since the magnetic field 
couples only to charged molecules and grains, which in 
turn determine the evolution of the magnetic field and 
thus the dynamics of the contraction. In a companion 
publication (Tassis et al. 2011, in preparation, hereafter 



Paper II), we discuss in more detail the interplay between 
chemistry and dynamics in magnetic models. 

Our paper is organized as follows. In fJ5] we discuss 
our dynamical models (hydrodynamic and MHD) of core 
evolution. The chemical network we implement is dis- 
cussed in *j3j Details on the parameter study we have 
conducted are given in [2J We present our results in iJ5] 
and we discuss our findings and conclusions in fj6] 

2. DYNAMICAL MODELS OF MOLECULAR CLOUD CORES 

2.1. Hydrodynamic Core Models 

We model the contraction of a non-magnetic, self- 
gravitating, isothermal prestellar core in spherical sym- 
metry. The core is assumed to be embedded within a 
uniform-density isothermal cloud. We solve the equa- 
tions of hydrodynamics in spherical symmetry: 
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where v r is the radial velocity, p is the total density of 
the gas, and pj is the density of each individual chemi- 
cal species followed (listed in Table [T]). We solve a den- 
sity continuity equation for each of the chemical species 
in the model (equation [TbJ) . In these equations, Sj are 
source/sink terms describing mass interchange between 
the species due to chemical reactions. 

The initial state is one of uniform density equal to 
10 3 cm -3 and zero velocity. The radius of the core is 
0.4 pc, making it significantly thermally supercritical. 
At the outer boundary, we require constant pressure and 
zero velocity. This way, the core experiences no influx of 



2.2. MagnetoHydroDynamic Core Models 

We follow Basu & Mouschovias (1994) in modeling 
the formation and contraction of prestellar cores in self- 
gravitating, magnetic, rotating, isothermal model clouds 
in the presence of ambipolar diffusion and magnetic 
braking. In the presence of an ordered magnetic field 
the cloud relaxes rapidly along field lines until thermal- 
pressure forces balance gravity, before any significant 
contraction perpendicular to the field lines takes place 
(Fiedler & Mouschovias 1993; Desch & Mouschovias 
2001; Kunz & Mouschovias 2010). Force balance along 
magnetic field lines is maintained even when the con- 
traction becomes dynamic in the perpendicular direction. 
Consequently, the thin-disk approximation is adopted for 
the model cloud cores (Ciolek & Mouschovias 1993). The 
model cloud is axisymmetric with the rotation axis coin- 
ciding with the axis of symmetry, which is taken to be the 
z-axis of a cylindrical polar coordinate system (r,(f>,z). 
The upper and lower surfaces of the disk are located at 
z = ±Z(r,t). The cloud is embedded in an "external" 
medium of constant thermal pressure P cx t, and density 
p ox t, and is threaded by a magnetic field which, far from 



the cloud, becomes uniform and equal to a "reference" 
value B rc {. 

The thin disk evolution equations, obtained by inte- 
grating the two-fluid (plasma and neutrals) MHD equa- 
tions over z, are 
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where er n and p n are the column and volume densities 
of neutrals respectively, Vi ir and v n ^ r are the radial ion 
and neutral velocities respectively, L n = a n v u>( f,r is the 
angular momentum per unit area, and $ = dr'r 1 ' B z cq 
is the magnetic flux. Equation (I2al) is the mass conti- 
nuity equation of the neutrals. Equation (|2c[) gives the 



4 



<f> component of the neutral force equation, while the r- 
component is given by Equation (|2b|) . Equation ((2cJ) 
states that the thermal pressure at the equatorial plane 
in each magnetic flux tube balances the external and 
gravitational pressures. The introduction of the effective 
sound speed, C Q g, accounts for the contribution of the 
radial force exerted by the external pressure P ex t on the 
upper and lower surfaces of the flaring disk, since the lat- 
ter are not horizontal. Furthermore, in the expression for 
the total radial magnetic force (per unit area), i^M.r, the 
first two terms are the magnetic tension force and the 
magnetic pressure force within the cloud, respectively. 
The remaining terms are corrections due to the variation 
of t he half-thickness of the cloud, Z, with r. Equation 
(|2d[) is Faraday's law in the two-fluid approximation with 
the magnetic flux frozen in the ion fluid. The toroidal 
component of the magnetic field (-B^z) is determined 
assuming a steady state condition on the transport of 
angular momentum by torsional Alfven waves in the ex- 
ternal medium and given by equation (j2"H) . The radial 
component of the magnetic field at the top surface of 
the thin disk (-B r ,z) is given by equation ( [2j| >, obtained 
by assuming a force-free external medium, whereby the 
half-thickness of the disk (Z) is obtained from equation 
(|2f|) after p n (r,t) is calculated from equation (|2c|) . The 
r-component of the gravitational field (|2kj) is calculated 
from Poisson's equation for a thin disk. Jo is the zeroth- 
order Bessel function of the first kind, K is the com- 
plete elliptic integral of the first kind, r< = min[r, r'], 
and r> = max[r, r']. Equation (|2mj) is the mass con- 
tinuity equation for neutral gas phase chemical species 
advected with the neutral bulk velocity u n ,r- Equation 
(|2n|) is the mass continuity equation for the gas phase 
ions and all the grain mantle species (since dust grains 
are overwhelmingly negatively charged in the densities 
we are studying) advected with the velocity of the ions 
v- ur . Sjo and Sj+ represent the source/sink terms due to 
chemical reactions. Note that the ions are advected with 
a different velocity than the neutrals, so two equations 
are required in the magnetic model to replace Eq. (jlb|) 
of the non-magnetic model. 

Equation (|2o]) gives the velocity of the ions and is ob- 
tained from the algebraic solution of the force equation 
for the ion fluid where the acceleration terms have been 
neglected because of the small inertia of these species. 
The mean (momentum exchange) collision time of a neu- 
tral particle with ions in equation (|2o[) is 

< r ni >= 1.56— 2 S 3) 

Pi < aw > in 

where m p is the proton mass, p\ = Pi+ 1S the total 
ion density, the mean molecular weight of the ions is 
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7ij + is the number volume density for each one of the 
ions considered, and < aw > ln is the average collisional 
rate between ions and hydrogen molecules. The value 
of < aw >i n varies slightly with the mass of the ion 
(McDaniel & Mason 1973); we use the canonical value 
< aw >in= 1-69 x 10 -9 cm 3 s _1 (Ciolek & Mouschovias 
1993). 

Each magnetic model is initiated at a reference state 
as in Basu & Mouschovias (1994). The reference state 



has a column density profile of the form 
a n ,ref{r) 



[l + (r/Z rcf ) 2 ] 3/2 ' 
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which implies an almost uniform column density for 
r <C Z re f, which is where the core forms. The high-r tail 
has little effect on the calculations and is there for nu- 
merical convenience, allowing the calculation of the disk 
gravitational potential to converge. The profile of the ref- 
erence state angular velocity is obtained by using a linear 
dependence of specific angular momentum on mass, cor- 
responding to the angular momentum distribution of a 
uniformly rotating disk, yielding 
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The reference state is first allowed to relax to an equilib- 
rium configuration, while ambipolar diffusion and chem- 
istry are turned off. It is from this equilibrium state that 
the subsequent evolution is initiated. For this reason, 
the initial central density to which each model relaxes is 
not always the same but it depends on the initial mass 
to flux ratio (see Fig. [1]). 

The free parameters in the magnetic models are the 
central density of the initial reference state (which 
we take to be 10 3 cm~ 3 ); the background angular 
velocity fit = 10 _15 rads _1 , approximately corre- 
sponding to the rate of galactic rotation in the so- 
lar neighborhood; the reference state angular velocity 
(O re f = 0.48 kms^ 1 pc^ 1 ); the external pressure P cxt = 
0.l7rGcr? ref /2 corresponding to 10% of the central "grav- 
itational" pressure; the characteristic size of the col- 
umn density profile l le f = 1 pc; the external density 
Pcxt = 10cm -3 ; and the mass to magnetic flux ratio, 
which is varied (see discussion of parameter study be- 
low). 

With the exception of the magnetic mass to flux ratio, 
the dynamical evolution of the collapsing core is not sen- 
sitive to any of the free parameters (Basu & Mouschovias 
1995). 

At the outer boundary of the disk we require continuity 
of thermal as well as magnetic pressures (see Tassis & 
Mouschovias 2005, §5). No influx of mass is allowed at 
the outer boundary, so the total mass of the cloud is 
conserved. 

3. CHEMICAL NETWORK 

We use a subset of the UMIST chemical reaction net- 
work (Woodall et al. 2007), which includes data on 420 
species and 13 elements followed in gas-phase reactions. 
Each chemical species is treated as a separate fluid in our 
formalism. The source and sink terms on the right-hand 
side of Eq. (|lb|) and in Eqs. (|2m[) and (|2n[) describe the 
interactions between chemical species. They are given 
by: 

c _ d Pj 

~^ _ 

(6) 

where otuj is the rate coefficient for the production of 
species j due to the interaction of species k and I, j3j}~ 
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TABLE 1 
Chemical Species Considered 



gas phase species 

H+ H W 2 He He+ C+ C CH CH+ CH+ CH~ 2 N N+ CH^ 

NH+ CH^ NH NH^ O ch 4 cnj 0+ NH 2 CH^ OH OH+ NH^ nh 3 h 2 o nhJ 

H 2 0+ H3O+ C 2 C+ C 2 H+ C 2 H C 2 H+ C2H2 CN CN+ HCN+ C 2 H+ HCN HNC Si+ C 2 H+ 

H 2 NC+ Si N 2 CO+ HCNH+ CO N+ HCO N 2 H+ HCO+ H 2 CO H 2 CO+ NO NO+ H3CO+ CH3OH 

Q 2 O^ CH3OH+ C+ C 3 H+ C 2 N+ CNC+ C 3 H+ CH 3 CN C3H2 CQ 2 CO+ HCO+ HC 3 N 

grain mantle species 

H C CO H 2 CO Si C~ 2 Ol CH OH NO CH~ 2 H^O CO~ 2 CH~ 3 CH~ 4 HNC 

C 2 H 2 HC 3 N N 2 CN NH HCN C 2 H NH 3 CH3CN CH3OH NH 2 N O H 2 HCO C 3 H 2 

CH2OH 



TABLE 2 
Initial Abundances 

"h 1.00 x icr 2 

H 2 4.95 x lO" 1 

He 1.40 x lCT 1 

C+ 7.30 x 10" s 

N 2.14 x 10" 6 

O 1.76 x 10~ 4 

Si 2.00 x 10~ 8 



is the rate of production of species j from species k by 
photoreactions, and 8ij is the Kronecker delta function. 

In our models the source terms are implemented in an 
operator-split fashion, and the change of abundances due 
to the source terms is updated after the advection terms. 
The species we consider are listed in Table 1, and they 
include 78 gas-phase atomic and molecular species, and 
33 grain ice mantle species, that are composed of the 
elements H, He, C, N, O, and Si as a representative low 
ionization potential metal. The chemical network follows 
1553 reactions. 

All the models studied in this work use the isothermal 
approximation (small deviations from isothermality have 
been shown to not affect the dynamical evolution of cores 
significantly, e.g. Keto & Caselli 2010). We adopt as a 
fiducial temperature T = 10K, consistent with line in- 
tensity ratios of different molecular transitions showing 
that starless cores are almost isothermal, with T ~ 10 
K, (Tafalla et al. 1998, 2002, 2004), although there may 
be small variations of the order of a few K (e.g., Evans 
et al. 2001). In our parameter studies, we test for the 
effect of small temperature variations on our chemical 
abundances. 

Since the prestellar cores we model are embedded in 
a molecular cloud we assume a visual extinction Ay = 
3 mag at the outer boundary, so that the results arc 
not much affected by photodissociation. For the inner 
layers of the core we obtain the visual extinction from 
the column density to the surface using the empirical 
conversion relation Nn 2 /Ay = 9.4 x 10 20 cm~ 2 mag -1 
(Bohlin et al. 1978). For our deeply embedded cores, 
no H 2 self-shielding is required; CO self-shielding on the 
other hand is accounted for based on the model of Lee 
et al. (1996). 

We assume a standard MRN distribution of grains 



(Mathis, Rumpl & Nordsieck 1977), assumed to be well- 
mixed with the gas. Since most of the grains are neg- 
atively charged at the densities of interest in prestellar 
cores, we also take into account the electrostatic attrac- 
tion between ions and grains through an enhanced rate 
of accretion for the ions (Rawlings et al. 1992). We 
use a fixed grain abundance (n gr ain/«ncutrai = 10~ 12 ) 
and charge (we take all grains to be singly negatively 
charged). Grain processes modeled include freeze-out, 
surface reactions, thermal desorption, and cosmic ray 
heating. We have used a sticking probability of 0.3. 

The discrete nature of grains means that chemistry on 
their surfaces should properly be modeled by stochastic 
methods (e.g.; Stantcheva, Shematovich, & Herbst 2002; 
Chang et al. 2007). However, those methods are compu- 
tationally very expensive, especially for a large chemical 
network and the large number of HD and MHD models 
we run. Instead, we use the rate equations method, with 
modifications as suggested by Vasyunin et al. (2009), 
which brings the results of this method into good agree- 
ment with Monte Carlo modeling for the range of tem- 
peratures we are concerned with. The grain reaction net- 
work is taken from Hascgawa & Herbst (1993), as are the 
rates for desorption by cosmic ray heating of grains (up- 
dated for new binding energies where necessary). Where 
available, we use binding energies determined in the lab 
(see Willacy 2007). Other binding energies are taken 
from Hasegawa & Herbst (1993). The CO2 chemistry on 
the grains is based on Garrod & Pauly (2011) 4 . 

The gas phase reactions are better constrained as they 
have been studied for a longer time, and there ap- 
pears to be convergence in the literature among differ- 
ent groups (e.g., Aikawa et al. 2003) with the exception 
of sulphur-bearing molecules (Wakelam, Herbst & Selsis 
2006), which, however, are not considered in our model. 
In our study we aim to assess the effect of the dynam- 
ics on chemical abundances, while the chemical network 
remains fixed. 

We use low-metal initial abundances 5 , given in Table 
2. The calculations are initiated with all elements other 
than molecular hydrogen in the form of atoms. Of course 
in nature there is some chemical evolution taking place 

4 with the activation barrier for the CO + O — > CO2 reaction 
on the grain ice increased from 290 K to 340 K (within the exper- 
imental range). 

5 with the exception of Si, the abundance of which is higher than 
in the low-metal model 
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during the formation of the parent cloud. A hrst attempt 
to calculate initial abundances for molecular clouds for 
one particular cloud formation scenario was presented by 
Hassel et al. (2010). However, these results are model- 
dependent and will differ for different cloud formation 
models. We will return to the issue of initial abundances 
in a future publication. 

4. PARAMETER STUDY 

For each class of models (non-magnetic or magnetic), 
we examine the dependence of the molecular abundances 
and their evolution on four parameters: the C/O ratio, 
the cosmic ray ionization rate, the temperature, and a 
parameter controlling the time available for chemical evo- 
lution. In non-magnetic models, the latter is a time in- 
terval prior to collapse during which chemistry evolves 
but the core is assumed to not evolve dynamically. In 
magnetic models, it is the initial mass-to-magnetic-flux 
ratio which controls the amount of magnetic support: a 
smaller initial mass to magnetic flux ratio implies a larger 
amount of magnetic support against gravity, and hence 
a slower dynamical evolution for the core. 

Our "reference" non-magnetic model has a C/O ratio 
of 0.4, a cosmic ray ionization rate of £ = 1.3 x 10~ 17 s~ l , 
a temperature T = 10 K and an initial collapse delay time 
of 1 Myr. Our "reference" magnetic model similarly has 
a C/O'ratio of 0.4, £ = 1.3 x lO^^s- 1 , T = 10 K, and a 
central mass-to-flux ratio equal to the critical value for 
collapse (Mouschovias & Spitzer 1976). 

For non-magnetic models, we examine two additional 
values of delay: zero, and 10 Myr. For magnetic mod- 
els, we examine two additional values of the initial mass 
to magnetic flux ratio: 1.3 times the critical value (a 
faster-evolving, magnetically supercritical model), and 
0.7 of the critical value (a slower, magnetically subcriti- 
cal model). The three values of the delay and the three 
values of the mass-to-flux ratio set our dynamical model 
variations. When these models are combined with refer- 
ence values for the C/O ratio, £, and T, we will refer to 
them as the six basic dynamical models. 

The carbon-to-oxygen ratio is varied from its reference 
value by keeping the abundance of C constant and chang- 
ing that of O (as in Terzieva & Herbst 1998). The two 
other values of C/O ratio examined are 1 and 1.2. We 
have varied each of the six basic dynamical models by 
changing its C/O ratio to one of the values above, and 
we have thus studied 12 models (6 magnetic and 6 non- 
magnetic) with the C/O ratio varied from its reference 
value. 

We have also considered different temperatures for 
each of the six basic dynamical models. We have varied 
the temperature by a factor of ~ 1.5 from its reference 
value of 10 K and examined models with T = 7 K and 
T = 15 K, and we have thus studied 12 models (6 mag- 
netic and 6 non-magnetic) with a temperature differing 
from its reference value. 

Finally, to test the effect of the cosmic ray ionization 
rate we have studied two additional magnetic and two 
non-magnetic models, which are, in each case, identical 
to the corresponding "reference" models, but with differ- 
ent cosmic-ray ionization rates: namely, a factor of four 
above (C = 5.2 x 10" 17 s" 1 ) and below « = 3.3 x 10~ 18 
s _1 ) its reference value (covering the range of observa- 
tional estimates e.g.: MaCall et al. 2003; Hezareh et al. 



2008). 

Thus, we have examined a total of 34 models (6 basic 
dynamical models, 12 with varied C/O ratios, 12 with 
varied T, and 4 with varied £). These 34 models are 
depicted in our figures in [J5] with different line types, 
thicknesses, colors, and symbols, as follows: 

• Line type (solid, dashed, dotted) and symbols (X, 
diamond, cross) depict the evolutionary timescale 
("reference", "slow", and "fast" respectively). 

• Color denotes the type of model (magnetic or non- 
magnetic), with bluish hues corresponding to non- 
magnetic models, and reddish hues corresponding 
to magnetic models. 

• The line thickness/symbol size (for red and blue 
lines and symbols) indicate the values of the cos- 
mic ray ionization rates: thin red/blue solid lines 
and small red/blue symbols correspond to the low 
value of £, while thick red/blue solid lines and large 
red/blue symbols to the high value of £, with all 
other parameters fixed at their "reference" values. 

• The effect of different C/O ratio values is shown as 
a shaded area (orange for magnetic runs, and cyan 
for non-magnetic runs). The edge of the shaded 
area in contact with the corresponding ("fast", 
"reference" , or "slow" , magnetic or non-magnetic) 
model line corresponds to C/O=0.4, while the edge 
of the shaded area furthest away from the line cor- 
responds to C/O =1.2 

• Finally, the effect of the temperature is shown 
with the colors brown for magnetic runs and pur- 
ple for non-magnetic runs, either as a shaded 
area, or as line thickness/symbol size. In the 
latter case, the thick brown/purple lines (or the 
large brown/purple symbols) correspond to T = 
15 K and the thin brown/purple lines (small 
brown/purple symbols) to T = 7 K. 

5. RESULTS 

5.1. Dynamics 

Figure [1] shows the evolution of the central density 
of contracting core models with time. Red lines corre- 
spond to magnetic runs, and blue lines to non-magnetic 
runs. The different magnetic runs correspond to dif- 
ferent values of the mass-to-magnetic-flux ratio: criti- 
cal for the solid lines/ x symbols; 0.7 times the critical 
value (magnetically subcritical, slower-evolving magnetic 
model) for the dashed lines/diamonds; and 1.3 times the 
critical value (magnetically supercritical, faster-evolving 
magnetic model) for the dotted lines/ + symbols. The 
different non-magnetic runs correspond to different as- 
sumed delay times, during which the chemistry oper- 
ates, before the core is allowed to evolve dynamically 
(collapse): 1 Myr delay for the solid lines/ x symbols; 10 
Myr delay (slowest-evolving non-magnetic model) for the 
dashed lines/diamonds; and no delay (fastest-evolving 
non-magnetic model) for the dotted lines/ x symbols. 
The symbols mark the times when full radial outputs 
are obtained and analyzed, and correspond to successive 
enhancements of the central density by an order of mag- 
nitude. 
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Fig. 1. — (a) Central number density, n c , as a function of 
time. Symbols indicate the times at which full radial snapshots 
are taken. Blue lines: non-magnetic models (dotted/crosses: no 
delay; solid/x: 1 Myr delay; dashed/diamonds: 10 Myr delay). 
Red lines: magnetic models (dotted/crosses: supercritical model; 
solid/x: critical model; dashed/diamonds: subcritical model). The 
thin and thick red solid lines correspond to cosmic ray ionization 
rate £ a factor of four above and below the "reference" value, and 
bracket the normal-thickness red solid line. The yellow shaded ar- 
eas show the effect of a varying C/O ratio, and correspond to a 
range of C/O values between 0.4 (reference value) and 1.2. The 
brown shaded areas show the effect of varying the temperature 
from 7K to 15K, the reference value being 10K. 



Different dynamical models differ with respect to the 
two fundamental ingredients affecting the chemical evo- 
lution of a core: density, and time spent at each density. 
Non-magnetic models with substantial time delays (rep- 
resenting the time it takes for some level of turbulent 
support to dissipate) can take as much time as substan- 
tially magnetically subcritical models to reach very high 
densities. However, once the collapse sets in (support 
dissipates), the evolution occurs essentially within one 
free-fall timescale, so that the bulk of the time delay is 
spent at the initial, low density (10 3 cm -3 ). In contrast, 
slow magnetic models allow the core to spend substan- 
tial time at higher densities, and as a result the chemical 
state of these models are expected to be different. 

The thin and thick red lines correspond to magneti- 
cally critical models but have ( a factor of four above 
and below the "reference" value, and they bracket the 
normal-thickness red solid line. Because the ambipolar 
diffusion timescale depends on the degree of ionization, 
which in turn depends on the cosmic-ray ionization rate, 
( affects the evolutionary timescale of the magnetic mod- 
els: a higher degree of ionization implies increased cou- 
pling between magnetic field and neutral fluid, and a 
resulting slower evolution, as magnetic forces are more 
effective in opposing gravity (see Paper II). 

In addition to the mass-to-flux ratio and the degree 
of ionization, the evolutionary timescale of the magnetic 
models is also affected by the C/O ratio. The effect of 
different values of the C/O ratio can be seen in Fig. [1] as 
the orange shaded area, which denotes a range of C/O 
values between 0.4 (reference) and 1.2. The reason for 



this dependence is that the C/O value affects the chemi- 
cal reaction network and consequently the degree of ion- 
ization at each density, resulting to a different evolution 
rate. In general, higher values of the C/O ratio result in 
slower evolution, and the effect is more pronounced for 
models with a low mass-to-flux ratio. 

Finally, the evolutionary timescale of the magnetic 
models also depends on the temperature. The effect of 
different T values is shown in Fig. Q] as the brown shaded 
area, which denotes a range of T values between 7 K 
and 15 K. The dependence is stronger for slower mod- 
els, and it enters through the dependence of the degree 
of ionization on temperature (see Paper II). The evolu- 
tionary timescale increases with increasing temperature, 
and the difference between 10 K and 15 K is significantly 
larger than the difference between 7 K and 10 K, al- 
though the fractional increase is comparable (about 50% 
in each case). 

Radial profiles of the density, velocity, and (edge-on for 
the magnetic runs) column density for the magnetic and 
non-magnetic reference models and for the different time 
snapshots corresponding to the x symbols in Fig. [1] are 
given in Fig. [2] Again the magnetic model is represented 
by red lines and the non-magnetic model by blue lines. 
Each snapshot corresponds to a time when the volume 
density is an order of magnitude higher than the pre- 
vious one. The time elapsed between the initialization 
of each simulation and each snapshot for the magnetic 
model is 3.53 Myr, 4.4 Myr, 4.6 Myr, 4.65 Myr. The cor- 
responding times for the non- magnetic model are 1.89 
Myr, 2.02 Myr, 2.05 Myr, 2.06 Myr. The first snapshot 
in both cases corresponds to the initial configuration of 
each model. 

Radial profiles of the volume density are shown in the 
left panel. These are power-law profiles with a flat inner 
region, the boundary of which is indicated with a circle 
for the magnetic runs and a square for the non-magnetic 
runs. Profiles of the infall velocity, in units of the sound 
speed, arc shown in the middle panel. Here the differ- 
ence between the magnetic and the non-magnetic models 
is most striking, as for the magnetic model the infall is 
always subsonic, while for the non-magnetic model the 
infall velocities are substantially supersonic throughout 
the evolution of the core and for most of its radial ex- 
tent. The right panel shows radial profiles of the column 
density, which have similar properties with the volume 
density profile. 

The masses of the different model cores are not ex- 
actly the same. This is because the mass of each core 
depends on the temperature in the case of non-magnetic 
cores, and the mass-to-flux ratio, temperature, C/O ra- 
tio, and cosmic-ray ionization rate in the case of mag- 
netic cores. In spite of these dependencies, we have se- 
lected the initial conditions of our models so that the core 
masses remain roughly comparable: the non-magnetic 
"reference" model core mass is 15 M Q and the magnetic 
"reference" model core mass is 20 M Q , while in general 
the core masses in each class of models do not differ from 
this value by more than a factor of ~ 2. 

A shared feature of the volume and column density ra- 
dial profiles with important observational implications is 
the size of the flat-profile inner region. It is well-known 
that prestellar cores have uniform column density inner 
regions (Ward- Thompson et al. 1994; Ward-Thompson et 




Fig. 2. — Radial profiles of: the number volume density, n (panel a); velocities of the neutrals in units of the isothermal sound speed 
(panel b); and (edge-on for the magnetic runs) column density (panel c). Red lines: magnetic runs; Blue lines: non-magnetic runs. Each 
curve corresponds to a snapshot in time when the volume density is an order of magnitude higher than the previous one. Time elapsed per 
snapshot in magnetic runs: 3.53 Myr, 4.4 Myr, 4.6 Myr, 4.65 Myr; in non-magnetic runs: 1.88 Myr, 2.01 Myr, 2.05 Myr, 2.06 Myr. 



al. 1999; Bacmann et al. 2000; Shirley et al. 2000; Schnee 
et al. 2010). The angular size of the column density pro- 
file fiat inner region, which is observationally measurable 
in well resolved cores, combined with an estimate for the 
distance can yield the physical size of the column density 
profile flat inner region, which is of the same size as the 
volume density profile flat inner region, as can be directly 
seen by comparing the left and the right panels of Fig. [51 
The latter can be used to estimate the central volume 
density of a collapsing core (Tassis & Yorke 2011). 

5.2. N on- equilibrium Chemistry, Production and 
Destruction of Commonly Studied Molecules 

The need for the inclusion of non-equilibrium chem- 
istry in our dynamical models is explicitly demonstrated 
in Fig. [3] Here we plot the ratio P/D between the total 
production and destruction rates of various molecules, 
as a function of density, at the center of our reference 
core models. The left panel corresponds to the mag- 
netic reference model, and the right panel to the non- 
magnetic reference model. The total production and 
destruction rates are defined as including all chemical 
processes leading to the production or destruction of a 
specific molecule respectively. Each line corresponds to a 
different molecule. Note that the y-axis in the two plots 
is logarithmic and extends over ten orders of magnitude. 
A large number of molecules are very much out of equi- 
librium {P/D is very far from 1), and they remain so 
throughout each simulation. In the magnetic run, where 
more time is available for the system to approach chem- 
ical equilibrium even at high densities, the majority of 
molecules eventually move closer to P/D = 1 ; in the non- 
magnetic model the reverse trend is observed, and as the 
density increases and the evolution timescale decreases, 
most molecules deviate further away from P/D = 1. 
However, even for the magnetic model the spread about 
P/D = 1 is substantial, as demonstrated by the inset on 
the left panel, which zooms in around central densities 
of 10 4 — 10 6 particles per cm 3 and P/D within 50% of 
unity displayed on a linear scale. In the inset, the values 
of P/ D exhibit a spread of more than 20% around unity. 

The effect of non-equilibrium chemistry becomes obvi- 



ous when we examine the dependence of central molec- 
ular abundances on central density for different dynam- 
ical models. If the chemistry was in equilibrium, then 
the central volume density would uniquely determine the 
central molecular abundances. This, however, is not the 
case, as we can see in Fig. [4j where we examine the cen- 
tral abundances, production, and destruction rates due 
to various reactions for one carbon-bearing (CO) and one 
nitrogen-bearing (NH3) molecule. The left column cor- 
responds to the non-magnetic reference model (1 Myr 
delay), and the right column corresponds to the mag- 
netic reference model (critical mass-to- flux ratio). The 
top row shows the evolution with central density (bottom 
axis) and with time (top axis) of the central abundances 
of CO (blue lines) and NH3 (purple lines). Solid lines cor- 
respond to gas-phase abundances and dashed lines cor- 
respond to abundances on grain ice mantles. The middle 
and bottom rows show the most important production 
(black lines) and destruction (red lines) reactions for the 
two molecules. The importance of each reaction is quan- 
tified by the percentage by which the specific reaction 
contributes to the total production/destruction rate of 
each molecule (y-axis). 

Comparing the top left and right panels, we can im- 
mediately see that the central abundances of the two 
molecules at the same central density are very different 
for different dynamical models, a result of the chemistry 
being out of equilibrium. Because each model spends 
different amounts of time at each density (see Fig. [IJ, 
non-equilibrium chemistry results to significantly differ- 
ent central abundances. Note that since the magnetic 
model evolves more slowly, the gas phase abundances of 
both molecules in the magnetic core model decrease more 
steeply with density, as the molecules freeze out on the 
grain mantles; correspondingly, there is a faster increase 
in the ice mantle abundance of each molecule. 

It is evident from the middle and bottom panels that 
any specific reaction reaches its maximum importance at 
different densities for different dynamical models. Freeze 
out onto grain mantles (solid red line) , which becomes in- 
creasingly important at higher densities, is a prominent 
example. Furthermore, for each molecule different reac- 
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Fig. 3. — Ratio of total (including all possible channels) production and destruction rates as a function of central density, for the magnetic 
(left panel) and the non-magnetic (right panel) reference models. The inset in the left panel shows a zoom-in of the density range between 
10 4 and 10 6 cm , with a linear y axis between P/D = 0.5 and 1.5. 
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tions dominate the production and destruction channels 
for different dynamical models. 

5.3. Evolution of central abundances 

Figures [5] and [6] show the evolution with central density 
of the central abundances of selected molecules in the gas 
phase. Figure [5] shows the effect of the dynamics, C/O 
ratio, and £ on the gas-phase abundances, while Fig. [6] 
shows the effect of the temperature. 

In Fig. [5J reddish hues correspond to magnetic mod- 
els and bluish hues correspond to non-magnetic mod- 
els. Different line types correspond to different dynami- 
cal models, with dashed lines representing the "slowest" 
models (subcritical magnetic model and 10 Myr delay 
non-magnetic model); solid lines representing the "in- 
termediate" models (critical magnetic model and 1 Myr 
delay non- magnetic model); and dotted lines represent- 
ing the "fast" models (supercritical magnetic model and 
no-delay non-magnetic model) . The shaded bands (cyan 
for non-magnetic models, orange for magnetic models) 
show the effect of the C/O ratio on the abundances (the 
edge of the band adjacent to a model line corresponds 
to C/O=0.4 and the edge furthest away from the line 
to C/0=1.2). The line thickness is representative of 
the value of the cosmic ray ionization rate, with normal 
thickness lines corresponding to £ = 1.3 x 10~ 17 s _1 , and 
thick (thin) lines corresponding to values of £ a factor of 
four above (below) that value. 

A general feature of the central gas-phase abundance 
evolution curves of Fig. [5] is their dependence on the 
dynamical evolution model. Magnetic models, which 
evolve more slowly between different central densities 
and spend more time at each value of the density, have 
an abundance peak at low densities, and generally ex- 
hibit a faster drop of the gas-phase abundance with in- 
creasing central density, due to efficient freeze-out onto 
grain ice mantles. As a result, the abundance evolution 
curves appear convex. In contrast, non-magnetic models 
evolve faster, chemistry has less time in which to oper- 
ate at each density, the abundances peak and drop at 
higher densities, and the non-magnetic abundance evo- 
lution curves are generally concave. The quantitative ex- 
planation of this behavior can be given by comparing the 
dependence of the freeze-out timescale and the collapse 
timescale on density for different models. The steep drop 
in molecular abundances (especially in carbon-bearing 
molecules) occurs at the density where the freeze-out 
timescale becomes shorter than the collapse timescale 
(Walmsley 1991). In magnetic models this happens at 
much lower densities (~ 10 4 cm -3 ) than in non- magnetic 
models (~ 10 6 cm~ 3 ). 

The results of Fig. [5] make it obvious that different col- 
lapse delay times at the initial phases of non-magnetic 
core evolution are not a good proxy for the effect of mag- 
netic fields. Although the slowest non-magnetic model 
curves (dashed curves, 10 Myr delay) are generally closer 
to the magnetic model curves, the difference is still signif- 
icant even in the qualitative behavior of the abundances 
(concave vs convex). The reason for this is that allow- 
ing more time for chemical evolution at ~ 10 3 cm -3 does 
not have the same effect, from a chemistry perspective, 
as spending a similar amount of time at higher densities, 
which is what happens in magnetic models. 

The behavior of certain molecules, such as H2CO, H^, 



NH, and CH3, is quite different for the non- magnetic 
models. While the slowest non-magnetic and all the mag- 
netic models have a qualitative similar behavior that re- 
sembles that of other molecules (the abundance builds 
up with increasing density, peaks at ~ 10 4 cm~ 3 , and 
then decreases at even higher density), the faster non- 
magnetic models, instead of a peak, exhibit a local min- 
imum at these densities. 

The molecules with the greatest sensitivity to the C/O 
ratio are, unsurprisingly, the oxygen-bearing molecules, 
the abundance of which is generally lower for higher C/O 
ratios. Molecules with more than one oxygen atoms are 
affected more than molecules with a single oxygen atom 
(see, e.g., CO versus CO2). Other molecules sensitive to 
the C/O ratio include CH, CH+, CH 3 , and NH. 

As expected, H^" responds directly to the change in 
the cosmic-ray ionization rate, as it is a direct ioniza- 
tion product (see Paper II for more details). Molecules 
that are primarily produced through reactions involv- 
ing direct ionization products, such as and Hc + 
(e.g., CH+) or free electrons, (e.g., NH 2 , NH, and CH 3 ) 
are similarly affected, even in the non-magnetic models. 
The magnetic models have an additional response to £ 
through the modification of their evolutionary timescale 
(see Fig. [IJ. For this reason, the change in molecular 
abundances, especially in the case of the largest- £ mag- 
netic model (which is also by far the slowest-evolving 
model) can be dramatic and feature very significant de- 
pletion in relatively low densities. The biggest difference 
is seen in molecules which, for moderate evolutionary 
timescales, show relatively modest depletion (e.g. N- 
bearing molecules). 

Figure [5] shows the effect of the temperature on the 
evolution with density of the central abundances of a 
few representative molecules. The blue and red lines 
correspond to the reference non-magnetic and magnetic 
models (the same ones shown with the normal-thickness 
blue and red lines in Fig. [5]). High-temperature (T = 15 
K) models are shown with the thick brown and pur- 
ple lines, for magnetic and non-magnetic models respec- 
tively, while low-temperature (T = 7 K) models are 
shown with the thin brown and purple lines; all other 
parameters are kept at their reference values. 

There are some significant differences between the 
molecular abundances in the 15K models compared to 
those in the 10K models, especially for nitrogen-bearing 
compounds. This is a result of two effects: (1) the reac- 
tion 

N+ + H 2 — ► NH+ + H (7) 

has a small activation barrier of 85 K. At 15 K the tem- 
perature is high enough for this barrier to be overcome 
and hence this reaction is much more efficient than at 
10K (see also Wakelam et al. 2010). This reaction is 
the start of the formation of nitrogen-bearing molecules, 
and an increase in its efficiency results in higher abun- 
dances of these molecules. (2) At 15 K the grains are 
warm enough for the desorption of N 2 ice to proceed more 
rapidly than at 10K (the rate of desorption is approxi- 
mately 10 10 times greater at 15 K compared to 10K). 

The way in which these processes act to produce the 
model results depends on the conditions in the models. 
As can be seen in Figure [51 increasing the temperature 
to 15K results in an increase in the abundance of ni- 
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Fig. 5. — Central abundances for representative molecules as a function of central density, for different dynamical models, different C/O 
ratios, and different cosmic ray ionization rates. Reddish hues correspond to magnetic models (dashed: magnetically subcritical; solid: 
magnetically critical; dotted: magnetically supercritical), and bluish hues to non-magnetic models (dashed: 10 Myr delay; solid: 1 Myr 
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"reference" value. The orange/cyan shaded areas show the effect of a varying C/O ratio in magnetic/non-magnetic models and correspond 
to a range of C/O values between 0.4 (reference value) and 1.2. 



trogen bearing molecules for the non-magnetic models 
(bold purple curves). The opposite effect is seen for the 
magnetic models (bold brown curves) . The difference be- 
tween the two types of models is due to the timescales 
over which they evolve. The non-magnetic models in Fig. 
[5] evolve quickest and reach their final density within 2 
Myr, whereas the magnetic models take between 5.3 and 
5.5 Myr (Fig. [IJ. Taking the non-magnetic models first: 
there is an increase in gaseous N2 caused by its faster 
desorption from the grains (see Fig. [5]). This means that 
the N2 abundance in the gas is higher at earlier times. 
The increase in N2 increases the destruction rate of H3 



as they react forming N2H + , and results in the sharp de- 
crease in H^" abundance seen in Fig.|6l Most of the N2H+ 
will be recycled back into N2 by recombination with elec- 
trons, but a small amount (~ 5%; Molek et al. 2009) will 
form NH. This can go on to form other gaseous molecules, 
or it can freeze out onto the grains, where it is quickly hy- 
drogenated to form ammonia ice. Ammonia has a much 
higher binding energy than N2 and once formed remains 
on the grains at these temperatures. The effect of the 
enhanced rate for reaction ([7]) can also be seen in the 
increased abundance of NH in the gas phase in the 15 K 
model. 
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Fig. 6. — Effect of temperature on central abundances of representative molecules. Blue/red lines: reference non-magnetic/magnetic 
models. Thick/thin brown lines: T = 15 K/ T = 7 K magnetic models; all other parameters kept at their reference values. Thick/thin 
purple lines: T = 15 K/ T = 7 K non-magnetic models; all other parameters kept at their reference values. 




The magnetic models are indicated by red/brown lines 
in Fig. |51 where the bold brown line is the 15 K model. 
The increased temperature here has exactly the oppo- 
site effect on nitrogen molecular abundances compared 
to the non-magnetic models. In this case the abundance 
of N2H + falls at 15 K, as do the abundances of most 
other nitrogen-bearing species. The abundance of H^" is 
roughly the same for all temperature models. The reason 
for this variance is the differing dynamical timescales of 
the models considered. The magnetic models remain at 



a relatively low density for > 2 Myr, before the onset 
of rapid collapse. During this prolonged period of time 
the chemistry evolves quite differently depending on the 
temperature. This is illustrated in Fig. [7] which shows 
the evolution of the N2H + abundance with time for two 
static models with a density of 2 x 10 4 cm~ 3 and with 
temperatures of 10 K (black) and 15 K (red). At 15 
K the nitrogen chemistry evolves faster, with X(N 2 H + ) 
reaching a peak value of > 10~ 13 cm~ 3 at 1 Myr, com- 
pared to 9 x 10~ 14 cm -3 at almost 2 Myr for the 10 K 
model. In the magnetic models, rapid collapse does not 
begin until after ~4 Myr. From Fig. [7] we can see that 
the N2H -1 " abundance is falling rapidly at this time and 
is more than an order of magnitude higher at 10 K than 
at 15 K. At 15 K much more nitrogen has been incor- 
porated into NH3 ices by this stage in magnetic models 
(see Fig. [8]) and less is available in the gas phase. Hence 
the early evolution of the collapse is quite different and 
results in the different abundances seen in Fig. [SJ 

Figure [5] shows the behavior of central abundances of 
various molecules in grain ice mantles, as a function of 
central density. Colors and lines are as in Fig. [SJ how- 
ever here the effect of temperature is also overplotted 
with a brown shaded area for the magnetic models, and 
a purple shaded area for the non-magnetic models. In all 
models, the composition of ice mantles is dominated by 
H2O, in accordance with observations of absorption fea- 
tures in the spectra of background stars viewed through 
dense clouds taken in IR (e.g., Whittet 2003 and ref- 
erences therein). In observed dark, quiescent regions 
(where cores are most likely to be prestellar and the ef- 
fect of the protostellar UV field is minimal), CO and 
CO2 are also observed to be important constituents of 
ice mantles: the contributions relative to the abundance 
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Fig. 8. — As in Fig. [5] for grain mantle ices. Reddish hues correspond to magnetic models (dashed: magnetically subcritical; solid: 
magnetically critical; dotted: magnetically supercritical), and bluish hues to non-magnetic models (dashed: 10 Myr delay; solid: 1 Myr 
delay; dotted: no delay). The thin and thick solid lines correspond to cosmic ray ionization rate f a factor of four above and below the 
"reference" value. The orange/cyan shaded areas show the effect of a varying C/O ratio for magnetic/non-magnetic models and correspond 
to a range of C/O values between 0.4 (reference value) and 1.2. The brown/purple shaded areas show the effect of varying temperature for 
magnctic/non-magnetic models and correspond to a range of T values between 7 K and 15 K (reference value is 10 K). 



of water ice is at the level of 10-40% for CO (Chiar et 
al. 1995; Gibb ct al. 2004; Chiar et al. 2011; Oberg et al. 
2011) and 20-30% for C0 2 (Whittet et al. 2009; Gibb et 
al. 2004; Oberg et al. 2011). Our models are in general 
agreement with these observations. There are, however, 
a few exceptions. The CO ice abundance in the 10 Myr 
delay time non-magnetic model shows a great sensitivity 
to temperature, with the 15 K model abundance being 
much lower (more than 3 orders of magnitude) than the 
rest of the models and observations. Similarly, the CO2 
abundance is decreased by much compared to observa- 
tions for the 7 K versions of the 10 Myr delay and 1 
Myr delay non-magnetic models, and all of the magnetic 
models. In observed clouds CH3OH is also found at the 
few % level relative to the water ice abundance (Boogert 
et al. 2011; Oberg et al. 2011); in our models the CH 3 OH 
abundance is at the order of < 1%. 

Nitrogen-bearing ices are dominated in our models 
by NH3 (< 10% of the water ice abundance), whereas 
CN ice is very underabundant, consistent with observa- 
tions (Gibb et al. 2004, Chiar et al. 2011). NH 3 ice is 
shown to be a good chronometer, with its abundance 
almost monotonically increasing with the evolutionary 



timescale. HCN is systematically more abundant than 
HNC. N 2 ice is very sensitive to temperature as discussed 
above, and its abundance decreases with increasing tem- 
perature. 

Molecular oxygen ice abundance has a sensitive depen- 
dence on the elemental C/O ratio (cyan and yellow bands 
for magnetic and non-magnetic models respectively) , and 
its abundance decreases as oxygen elemental abundance 
decreases, as expected. It is this dependence that dom- 
inates the range in possible values for the O2 ice abun- 
dance. 

The evolutionary trends with central density followed 
in our models by various mantle ice abundances are gen- 
erally similar: the abundance builds up fast as a result of 
efficient molecule freeze-out onto grains and high enough 
(Ay > 3) extinction, and saturates at higher values of 
the central density. 

5.4. Radial Abundance Profiles 

Figure [9] shows radial profiles of gas-phase abundances 
of representative molecules, taken when the core has a 
central density of 10 4 cm -3 (upper row) and 10 6 cm -3 
(lower row). Colors and lines are the same as in Fig. [8] 
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Fig. 9. — Radial profiles of gas- phase abundances of three commonly observed molecules, taken when the core has reached a central 
density of 10 4 cm -3 (upper row) and 10 6 cm -3 (lower row). Reddish hues correspond to magnetic models (dashed: magnetically subcritical; 
solid: magnetically critical; dotted: magnetically supercritical), and bluish hues to non-magnetic models (dashed: 10 Myr delay; solid: 1 
Myr delay; dotted: no delay). The thin and thick solid lines correspond to cosmic ray ionization rate £ a factor of four above and below the 
"reference" value. The orange/cyan shaded areas show the effect of a varying C/O ratio for magnetic/non-magnetic runs, and correspond 
to a range of C/O values between 0.4 (reference value) and 1.2. The brown/purple shaded areas show the effect of varying temperature for 
magnetic/non-magnetic models and correspond to a range of T values between 7 K and 15 K (reference value is 10 K). 



At high densities (lower row), both magnetic and non- 
magnetic models have experienced significant central de- 
pletion in all molecules (although generally to different 
extent). This is consistent with the qualitative picture 
of Fig. [5j at high densities the molecular abundances in 
all models have decreased significantly. At lower den- 
sities however (upper row), magnetic models exhibit a 
qualitatively different behavior from non-magnetic mod- 
els, at least for certain molecules. The difference is most 
striking in the case of CO, which for non- magnetic mod- 
els, is centrally-peaked at a central density of 10 4 cm~ 3 . 
In contrast, at the same density, magnetic models have 
already suffered significant depletion at the center. In 
the case of NH3 and CN, however, the radial profiles are 
quite similar (centrally peaked or flat) in both magnetic 
and non-magnetic models, even at low densities. As we 
will see in i )5.6[ this effect can be used to discriminate 
between non-magnetic and magnetic models. 

An important point to note is that the plots in Fig. [9] 
are in a logarithmic scale, and, as a result, most of the 
core mass is located in the larger scales depicted in the 
plots. For this reason, even in the cases where some or 
all models have undergone significant central depletion, 
when a core is not resolved the observed abundances sam- 
ple mostly the outer layers where little if any depletion 
has occurred. This results in degeneracies among models 
with otherwise significant differences in central depletion, 
when a core is unresolved. This point is prominently 



demonstrated in £15.51 

5.5. Total Abundances - Unresolved Cores 

Figure [TU] shows the total gas-phase abundances of se- 
lected molecules in the entire core (i.e. abundances mass- 
averaged over the entire extent of the "unresolved" core) , 
as a function of evolutionary stage, quantified by the 
core central density (x-axis). More advanced evolution- 
ary stages are to the right side of the horizontal axis. 
The effect of C/O ratio is shown with the shaded or- 
ange and cyan regions (for magnetic and non-magnetic 
models respectively), while the effect of temperature is 
shown with thin/thick brown and purple lines (thin for 
T = 7 K, thick for T = 15 K, brown for magnetic mod- 
els, purple for non- magnetic models). The effect of £ is 
only shown for the reference magnetic and non-magnetic 
models, and is represented by the thickness of the solid 
blue and red lines (thicker line corresponds to higher Q. 

A striking feature of the total gas-phase abundances 
is that there is little evolution with central density. The 
cause of this effect is the dominant contribution of the 
outer layers of the core (the lower-density regions) to 
the total core mass and their weight in determining the 
total core abundances as would be observed in unresolved 
cores. Since the abundances in the outer layers of the 
cores do not evolve much, the total core abundances also 
remain relatively stable with evolutionary stage. 

Because of the power-law abundance profiles, this fea- 
ture is generally preserved even if one averages only over 
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Fig. 10. — Abundances mass-averaged over the entire core for representative molecules at four different stages in the core's evolution 
(quantified by the core central density), for different dynamical models, different C/O ratios, and different cosmic ray ionization rates. 
Reddish hues correspond to magnetic models (dashed: magnetically subcritical; solid: magnetically critical; dotted: magnetically super- 
critical), and bluish hues to non-magnetic models (dashed: 10 Myr delay; solid: 1 Myr delay; dotted: no delay). The thin and thick solid 
lines correspond to cosmic ray ionization rate £ a factor of four above and below the "reference" value. The orange/cyan shaded areas 
show the effect of a varying C/O ratio for magnetic/non-magnetic models, and correspond to a range of C/O values between 0.4 (reference 
value) and 1.2. The brown/purple lines show the effect of varying temperature for magnetic/non-magnetic models. Thin brown/purple 
lines: T = 7 K; thick brown/purple lines: T = 15 K (reference value is 10 K). 



radii with densities above the critical density for exci- 
tation of a specific molecular line instead of over the 
entire core. Of course, the fact that a specific line 
is observed already places a lower limit on the cen- 
tral density of the core, equal to n cr it for the specific 
line. However, little evolution of the abundance with 
time is seen beyond that limit, unless n c rit is very high 
(> 10 6 cm~ 3 ). This can be seen in Fig. [TT] , where we 
have plotted abundances of selected molecules averaged 
over regions with densities higher than the critical den- 
sity n cr ;t for the following molecular lines: CO(J=l— >0), 
115 GHz, n > n CT it = 2.18 x 10 3 cm- 3 ; HCO+(J=1^0) 
89.2 GHz, n > n clit = 1.63 x 10 5 cnr 3 ; NH 3 (JK=11), 
23 GHz, n > 10 5 cm -3 [Although the critical density for 
this transition is much lower (~ 2 x 10 3 cm~ 3 ), this 
line requires much higher densities before the line is 
thermalized (Evans 1989)]; N 2 H+(J=1^0), 93.17 GHz, 
n > n clit = 2.25 x 10 5 c m - 3 ; HCN(JF=12^01), 88.631 
GHz, n > n clit = 2.18 x 10 6 c m - 3 ; and CN(N=l->-0), 
113.17 GHz, n > n clit = 1.3 x 10 6 cm" 3 . These criti- 



cal densities are for T — 10 K. All molecular data are 
taken from the Leiden Atomic and Molecular Database 
(Schoier et al. 2005). 6 

The large range of possible abundance values for differ- 
ent models at a fixed evolutionary stage (central density) 
demonstrates the large potential of chemical differentia- 
tion in prestellar cores. Even cores that appear identical 
in, e.g., their total column density profiles, can exhibit 
large variations in their molecular abundances as a result 
of a diversity in other local physical parameters (temper- 
ature, ionization rate, initial elemental abundance ratios, 
and dynamics). 

However, the large variation in molecular abundances 
among different models is neither monotonic nor system- 
atic, and the abundances resulting from different param- 
eters are not well-separated. Instead, we often see a non- 
linear dependence on certain parameters (a characteristic 
example is the ionization rate), as well as largely over- 

6 http:/ /www. strw.leidenuniv.nl/~moldata/ 
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Fig. 11. — As in Fig. 1101 but mass-averaged only over regions of the core with densities above the critical density n cr j t for excitation 
of commonly observed lines for each molecule: CO(J=l-H)), 115 GHz, n > n crit = 2.18 X 10 3 cm" 3 ; HCO+(J=l-s-0) 89.2 GHz, n > 
n crit = 1.63 X 10 5 cm- 3 ; NH 3 (J,K=1,1), 23 GHz, n > 10 5 cm- 3 , see text; N 2 H+(J=l->0), 93.17 GHz, n > n crit = 2.25 X 10 5 cm" 3 ; 
HCN(JF=12^01), 88.631 GHz, n > n crit = 2.18 X 10 6 cm" 3 ; CN(N=l->0), 113.17 GHz, n > n crit = 1.3 X 10 6 cm" 3 . 



lapping ranges, resulting in strong model degeneracies. 
We conclude that discrimination between models is not 
straight-forward, but instead requires the use of a variety 
of different observations and constraints, in combination 
with a statistical study accounting for as many different 
molecular abundances as possible. In the next sections, 
we will discuss certain quantities, such as depletion mea- 
sures and abundance ratios, which also hold promise as 
tools for the discrimination among dynamical models and 
physical parameters. 

We caution the reader that the abundances in Fig. [TU] 
are not to be compared directly with abundances derived 
from observations beyond an uncertainty of a factor of at 
least several. There are various reasons for possible dis- 
crepancies. First of all, the mass-averaged abundances 
presented here are not appropriate when the core is at 
least partly resolved; in this case, a more detailed study 
is needed, starting with the radial profile at a specific 
evolutionary snapshot for each model, and properly ac- 
counting for the beam size and the physical scale it cor- 
responds to at the distance of the core. Second, chemical 
reaction rates have appreciable uncertainties (Wakelam 
et al. 2010). We have not studied the effect of such un- 
certainties here, as our purpose was to investigate the 
effect of varying the core dynamics and other physical pa- 
rameters on molecular abundances. Finally, obtaining a 
molecular abundance from observations involves several 
approximations as well as assuming a value for the H2 
column density (which is not always well-constrained). 
Nevertheless, we have verified, for the following cores, 



that molecular abundances obtained from observations 
are generally within our model ranges, or only deviate 
from it within uncertainties: L1498 (Tafalla et al. 2004, 
2006); L1517B (Tafalla et al. 2004, 2006; Crapsi 2005); 
L1544 (Caselli et al. 1999; Jorgensen et al. 2004); L134N 
(Wakelam et al. 2006); TMC-l(CP) (Smith, Herbst & 
Chang 2004); L1689B (Jorgensen et al. 2004). For an 
improved comparison with observations, the appropri- 
ate methodology would be to start from our dynami- 
cal/chemical models, perform radiative transfer calcula- 
tions, and directly compare predicted line intensities with 
observed ones. We plan to return to such a comparison 
in a future study. 

5.6. Depletion Measures 

The qualitatively different behavior between the ra- 
dial profiles of magnetic and non-magnetic models seen 
in Fig. |9] can be used to discriminate between models. 
We define the central deviation A (central depletion or 
enhancement) of the molecular abundance as the ratio 
between the integrated abundance through the central 
30% 7 of the radial extent of the core over the integrated 
abundance over the entire core. 

Figure [12] shows the dependence of A on evolutionary 
stage (here quantified by the value of the central den- 
sity) for all models (colors and lines as in Fig. The 
horizontal dotted line marks the value A = 1. Values of 
A above the line indicate that the particular molecule is 

7 We have confirmed that the results are similar if this fraction 
is changed to 10% or 20%. 
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Fig. 12. — Central deviation of molecular abundance with respect to the average abundance over the entire core (A, see text), as a 
function of central density. The dotted line corresponds to A = 1 (no central enhancement or depiction). Reddish hues correspond to 
magnetic models (dashed: magnetically subcritical; solid: magnetically critical; dotted: magnetically supercritical), and bluish hues to 
non-magnetic models (dashed: 10 Myr delay; solid: 1 Myr delay; dotted: no delay). The thin and thick solid lines correspond to cosmic 
ray ionization rate f a factor of four above and below the "reference" value. The shaded areas show the effect of a varying C/O ratio, 
and correspond to a range of C/O values between 0.4 (reference value) and 1.2. The brown/purple shaded areas show the effect of varying 
temperature for magnetic/non-magnetic models and correspond to a range of T values between 7 K and 15 K (reference value is 10 K). 



centrally enhanced. Values of A close to the line indicate 
a flat abundance profile, while values of A below the line 
correspond to centrally depleted molecules. 

Despite the significant overlap between magnetic and 
non-magnetic models both in the radial profiles and 
the total abundances discussed in Figs. [9] and [10] re- 
spectively, for many molecules the magnetic and non- 
magnetic models not only are well separated, but also 
show little sensitivity to model parameters such as the 
evolutionary timescale, the temperature, the C/O ra- 
tio, and the cosmic-ray ionization rate. Molecules which 
show such good separation between magnetic and non- 
magnetic models include CO, C0 2 , HCO+, NO, OH, 
H2), and HC3N. In contrast, and CN show a par- 
ticularly degenerate behavior between different models. 

For non-magnetic models, CO-bearing molecules, and 
in particular CO, CO2 and HCO + , tend to show no de- 
pletion (A > 1) at early stages of the collapse, and only 
become depleted when the central density reaches 10 5 - 
10 6 cm~ 3 . In contrast, for magnetic models, the same 



molecules are already centrally depleted by central den- 
sities of ~ 10 4 cm~ 3 . 

NH3 on the other hand shows no sign of central deple- 
tion at central densities of ~ 10 4 cm~ 3 for both mag- 
netic and non-magnetic models. It only becomes de- 
pleted at later stages of the collapse (10 6 - 10 7 cm~ 3 for 
non-magnetic models, and 10 5 - 10 6 cm -3 for magnetic 
models). N2H + shows a similar behavior, however its be- 
havior is complicated by its sensitivity to temperature in 
the case of non-magnetic models (see also discussion in 
£ !5.3j) . Note also that in all N-bearing molecules, the high 
cosmic ray ionization rate magnetic model is always de- 
pleted, as a result of the very long evolutionary timescale 
(see also discussion of Fig|5]). This difference in behav- 
ior between CO and NH3 is consistent with observations 
(e.g., Caselli et al. 1999; Bacmann et al. 2003; Crapsi et 
al. 2007). 

5.7. Abundance ratios of molecules with discriminatory 

potential 
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Fig. 13. — Gas-phase abundance ratio between HCN and OH, 
mass-averaged over the entire core, which has good potential for 
discriminating between C/O ratio values. Red/blue symbols cor- 
respond to magnetic/non-magnetic models with T = 10 K. Di- 
amonds': "slow" models; X: "reference" models; crosses: "fast" 
models. The size of the red/blue symbols represents the value of 
f (larger symbols correspond to larger f). Brown/purple symbols 
correspond to magnetic/non-magnetic models with T = 7 K (small 
symbols) and T = 15 K (large symbols), and a "reference" value for 
f . The C/O ratio of each model is shown on the vertical axis, while 
the value of the abundance ratio is shown on the horizontal axis. 
Points correspond to a central H2 density of 10 6 cm -3 , however the 
core-average abundance ratios evolve little with density. 



Ratios between abundances of different molecules have 
been long considered as promising tools to overcome 
model uncertainties and degeneracies. However, most 
abundance ratios exhibit large variations and no consis- 
tent trend with any single model parameter when other 
parameters are varied. We have identified several gas- 
phase abundance ratios (averaged over the entire core) 
that are most affected by a single model parameter and 
show little variation as a result of the rest and thus show 
promising discriminatory potential. We plot these in 
Figs. [13] and [H 

Figure [13] shows the values of the abundance ratio be- 
tween HCN and OH, mass-averaged over the entire core, 
for our various models. This ratio can serve as an in- 
dicator of the C/O ratio. Different models are plotted 
with different symbols as discussed in the caption. The 
C/O ratio, shown in the vertical axis, exhibits a trend 
(increases with increasing value of each abundance ra- 
tio). Similarly, Figure [14] shows three abundance ratios 
that trace the evolutionary timescale of a core. These ra- 
tios are X NH3 /X C o, X NH2 /X C o, and X N h 3 /X H co+- All 
the ratio values are taken at the time when the central 
density in each model has reached 10 6 cm -3 . However, 
the averaged abundance over the entire core shows little 
time evolution, as shown in Fig. 1101 

Each abundance ratio examined in these two figures 
shows an appreciable scatter at a fixed value of the pa- 
rameter that it can help constrain (C/O ratio or evo- 
lutionary timescale, respectively). However, a general 
trend is obvious, and there are significant parts of the 
C/O or timescale parameter space that can be excluded 



with measurements of each ratio. 

6. DISCUSSION AND CONCLUSIONS 

We have combined dynamical and non-equilibrium 
chemical modeling of both magnetic and non-magnetic 
contracting molecular cloud cores, in order to predict 
molecular abundances as a function of core evolutionary 
stage, and identify observables with high promise in dis- 
criminating among core models. 

To this end, we have undertaken an extensive param- 
eter study, following a total of 34 different models, and 
varying: 

(a) the core evolutionary timescale, controlled by the col- 
lapse delay time, in the case of non- magnetic models, and 
the initial mass-to-magnetic-flux ratio, i.e. the amount 
of magnetic support against gravity, in the case of mag- 
netic models; 

(b) the elemental C/O ratio in the core; 

(c) the cosmic-ray ionization rate, which affects the 
chemistry in both magnetic and non-magnetic models, 
but also the evolutionary timescale in magnetic models, 
by affecting the extent of coupling between the magnetic 
field and the neutral fluid; and 

(d) the temperature. 

Our main results can be summarized as follows. We 
have explicitly demonstrated that many reactions in our 
chemical network are out of equilibrium during the dy- 
namical evolution of the core, necessitating the inclusion 
of non-equilibrium chemistry and the coupled modeling 
of chemistry and dynamics to properly follow the molec- 
ular abundances. 

For different dynamical models, the molecular abun- 
dances at the center of the core are not the same, even 
for identical values of the central number density; rather, 
the amount of time that a model core has spent at each 
density is what determines the actual molecular abun- 
dances at each evolutionary stage. Therefore it is not 
possible to properly account for magnetic effects in core 
chemistry by using a collapse retardation factor in non- 
magnetic models. Even if the overall final age of the 
core is the same, the fact that the non-magnetic core 
model has spent very different fractions of this age at 
various values of the central density than a magnetic 
core model has, results in significantly different molec- 
ular abundances. 

Abundances in magnetic models peak at much lower 
densities and decrease much faster at higher densities. 
The evolution of grain mantle ice abundances is qualita- 
tively similar in magnetic and non-magnetic models, with 
a fast initial increase due to freeze-out, and an abundance 
saturation at higher densities and later times. The value 
of that saturation however is different in different mod- 
els. We find that the abundance of NH3 ice scales roughly 
with the evolutionary timescale of the model, varying by 
a factor of 4 for our model range, which corresponds to 
an evolutionary timescale range of about 10 Myr. NH3 
ice is thus a potentially important chronometer. 

Radial profiles of gas-phase abundances reveal that the 
most significant difference between magnetic and non- 
magnetic models is at the central parts of the core, while 
most of the core extent tends to evolve much less and 
is much more similar among different models. However, 
as only a small fraction of the mass is concentrated at 
the central layers of the core, the central differences are 
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Fig. 14. — Gas-phase abundance ratios, mass-averaged over the entire core, with high potential for discrimination among evolutionary 
timescales for the core. Diamonds': "slow" models; X: "reference" models; crosses: "fast" models. Red/blue symbols correspond to 
magnetic/non-magnetic models with T = 10 K and C/O =0.4. The size of the red/blue symbols represents the value of f (larger symbols 
correspond to larger Q. Orange/cyan points correspond to magnetic/non-magnetic models with T = 10 K, C/0=1.2, and a "reference" 
value for f. Brown/purple symbols correspond to magnetic/non-magnetic models with T = 7 K (small symbols) and T = 15 K (large 
symbols), C/O=0.4, and a "reference" value for £. The time it takes the core to reach a central density of 10 e cm~ 3 is shown on the vertical 
axis, while the value of the abundance ratio is shown on the horizontal axis. Points correspond to a central H2 density of 10 6 cm -3 , however 
the core-average abundance ratios evolve little with density. 



smeared out in unresolved cores. As a result, the evolu- 
tion of the total, mass-averaged molecular abundances 
with central density in any particular model is much 
milder, if at all present. The situation is similar with 
the total abundances of grain mantle ices. 

Differences in the C/O ratio tend to affect mostly the 
abundance of the oxygen-bearing molecules, while the 
cosmic-ray ionization rate £ has its largest effect on the 
nitrogen-bearing molecules. In the case of magnetic mod- 
els, the direct effect of the cosmic-ray ionization rate on 
chemistry is compounded with its indirect effect on core 
evolution timescale, and as a result magnetic models are 
much more sensitive to changes in £. 

The many, sometimes competing, sometimes com- 
pounded, effects of different model parameters on molec- 
ular abundances result in severe degeneracies among dif- 
ferent models. For this reason, we have identified molec- 
ular abundance ratios which may help to constrain, even 
in unresolved cores, specific model parameters: the evo- 
lutionary timescale (ratios between NH3 and CO; NH2 
and CO; NH 3 and HCO+) and the C/O ratio (ratio be- 
tween HCN and OH). 

In addition, we have demonstrated that the cen- 
tral deviation A (central depletion or enhancement) of 
molecules can be another useful tool in model discrimi- 
nation. While for many molecules values of A are degen- 
erate among different models, there are several molecules 
for which there is a clear separation between magnetic 
and non-magnetic models. For these molecules the value 
of A even shows very little sensitivity to other model 
parameters. Measurements of the central depletion of 
such molecules therefore have an excellent potential for 
constraining dynamical models of core evolution. For 
this reason, it is very important to place adequate em- 
phasis into understanding differential depletion and grain 
chemistry effects as much as possible, and thus gain con- 
fidence in the chemistry part of these calculations, which 
still features significant uncertainties (eg., Wakelam et 
al. 2010). 

The dynamical models we have used in this study have 
low dimensionality: the non-magnetic models are spher- 



ically symmetric (one-dimensional), while the magnetic 
ones are axially symmetric, and integrated over the z- 
direction (1.5 dimensions). We note that as a result 
of the distinct geometries, the magnetic model does not 
converge to the non-magnetic model in the limit of zero 
magnetic field. A zero magnetic field disk model can- 
not be followed with the current numerical scheme, as 
the thin-disk approximation breaks down in the absence 
of a magnetic field. Although higher-dimensional calcu- 
lations are desirable (e.g., van Weeren et al. 2009), it 
is currently not possible to perform a parameter study 
as extensive as the one undertaken here using higher- 
dimensional models: a single model in this study (in- 
cluding non-equilibrium chemistry) requires more than 
800 CPU hours for the magnetic runs and more than 
200 CPU hours for the non-magnetic runs. We are how- 
ever in the process of implementing true 2D calculations 
for select cases presented here, to assess the effect of low 
dimensionality on our calculations. 

Because the model cores are not multi-dimensional, 
multiplicity or fragmentation is by definition neither ac- 
counted for nor allowed. This however is not a se- 
vere disadvantage in our case, as multiplicity at these 
early stages of collapse and for the low densities consid- 
ered here (up to 10 7 cm~ 3 ) is not observed (Schnee et 
al. 2010). 

Finally our work has not accounted for turbulent mix- 
ing; the flows in our model are laminar. Any turbu- 
lent mixing would result in making the radial profiles 
shallower (Xie et al. 1995). However, dense molecular 
cloud cores are generally not observed to have significant 
amounts of turbulence since their line widths are approx- 
imately sonic (Myers 1983; Barranco & Goodman 1998; 
Goodman et al. 1998; Kirk, Johnstone, & Tafalla 2007). 
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APPENDIX 
CODE TESTS 

As a test of the dynamical part of our code we have checked that, for the (more complex) magnetic model, our code 
returns the same results as Basu & Mouschovias (1994) when identical initial conditions are used. 

To confirm that the chemistry implementation in our code behaves as expected, we have checked that the elemental 
abundances are conserved in our simulations. Figure IA15I shows the fractional changes in the abundance of the most 
underabundant element Si (abundance of 10~ 8 , see Table [2]) as a fraction of central volume density, for magnetic 
(red lines) and non-magnetic (blue lines) of different mass-to-magnetic-flux ratios and initial delays, respectively. The 
fractional change always remains < 10~ 3 for all simulations. We thus confirm that the codes conserve elemental 
abundances better than one part in a thousand of the total abundance even for the most underabundant elements. 




io 4 io 5 io 6 io 7 
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Fig. A15. — Fractional change of total elemental abundance of Si (xsi) with respect to the initial one (xsi a ) versus central density in 
collapsing core models. Pure hydrodynamical model cores are shown with blue symbols and lines, while MHD core models are shown with 
red symbols and lines as in Fig. [T] 



In addition, we have performed a test on the advection of chemical species. Retaining identical initial and boundary 
conditions as in the main runs, in both magnetic and non-magnetic models, we allowed the chemistry to evolve for 
one million years before any dynamical evolution was allowed to take place. In this way, an abundance profile for the 
various molecules was developed. During the next step of the test, the dynamics were turned on, while the chemical 
reactions were turned off, so that molecules were advected, but their abundances could not be altered due to sources 
or sinks. Advection should not be changing the distribution of abundance with mass, and the abundances within a 
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Fig. A16. — Fractional change of the CO abundance, as a function of mass, at different times of the collapse, normalized to the initial 
time to. Left panel: non- magnetic reference model; Right panel: magnetic reference model. The inset shows the initial abundance profile. 
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single mass element should be conserved. The only expected change, as the contraction of the core proceeds, is the 
distribution of mass elements with scale. The results are presented in Fig. IA16I The left panel corresponds to the 
non-magnetic reference model and the right panel to the magnetic reference model. The plots show the relative error 
(normalized to the initial time) of the abundance of CO for different mass elements for different snapshots in time 
after the dynamics are turned on. The CO abundance profile is shown in the insets. We conclude that no spurious 
numerical alteration in molecule abundances is induced because of advection, as expected for a fully conservative 
advection scheme such as the one implemented here (Morton, Mouschovias, & Ciolek 1994). 
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